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ABSTRACT 

We study the nonlinear evolution of the magnetic Rayleigh-Taylor instability using three-dimensional MHD 
simulations. We consider the idealized case of two inviscid, perfectly conducting fluids of constant density sepa- 
rated by a contact discontinuity perpendicular to the effective gravity g, with a uniform magnetic field B parallel 
to the interface. Modes parallel to the field with wavelengths smaller than Ac = R - B/iph- pi)g are suppressed 
(where p/, and pi are the densities of the heavy and light fluids respectively), whereas modes perpendicular to 
B are unaffected. We study strong fields with varying between 0.01 and 0.36 of the horizontal extent of the 
computational domain. Even a weak field produces tension forces on small scales that are significant enough to 
reduce shear (as measured by the distribution of the amplitude of vorticity), which in turn reduces the mixing 
between fluids, and increases the rate at which bubbles and finger are displaced from the interface compared to 
the purely hydrodynamic case. For strong fields, the highly anisotropic nature of unstable modes produces ropes 
and filaments. However, at late time flow along field lines produces large scale bubbles. The kinetic and magnetic 
energies transverse to gravity remain in rough equipartition and increase as t'^ at early times. The growth deviates 
from this form once the magnetic energy in the vertical field becomes larger than the energy in the initial field. 
We comment on the implications of our results to Z-pinch experiments, and a variety of astrophysical systems. 



1. INTRODUCTION 

When a light fluid accelerates (or supports against gravity) 
a heavier fluid, the interface between the two is subject to the 
Rayleigh-Taylor instability (RTI). The idealized case of two in- 
compressible, inviscid, unmagnetized fluids of constant density 
separated by a contact discontinuity perpendicular to the ef- 
fective gravity g has been extensively studied through theory', 
experiment^-^ and numerical simulation^. In the linear regime, 
short wavelength perturbations of the interface grow fastest'. 
Once the perturbation amplitude is comparable to the wave- 
length, the interface can be characterized as rising "bubbles" 
of light fluid between descending "fingers" of heavy fluid. Sec- 
ondary Kelvin-Helmholtz instabilities are induced by the shear 
between the rising and descending columns. In the fully non- 
linear phase of a multimode spectrum of perturbations, smaller 
bubbles merge into larger structures, which then break up due 
to secondary instabilities, and a turbulent mixing layer results. 

From self-similarity arguments^, the time evolution of the 
height of the bubbles h above the initial interface is expected to 
be 



h = oAgt 



(1) 



where a is a dimensionless constant, and A is the Atwood num- 
ber 

(2) 

Ph + Pi 

iph and pi are the densities of the heavy and light fluids re- 
spectively). Recently, a comparison of the values of a mea- 
sured from laboratory experiments with the results of high res- 
olution, three-dimensional numerical simulations of multimode 
RTI with strong mode coupling, computed with a variety of nu- 
merical methods, has been presented by Dimonte et al^ (see 



also Ref. 4) as part of a validation effort for numerical meth- 
ods for hydrodynamics. Detailed analysis of the self-similar 
bubble dynamics, energy balance, and spectral properties of 
the resulting turbulent mixing layer demonstrated there is rea- 
sonable agreement between the simulations, theory and exper- 
iment, except in that the experimentally determined value of a 
is 0.057 ± 0.008, whereas most of the numerical simulations 
give a value for a which is about a factor of two smaller. It 
would appear this discrepancy is primarily due to mixing at 
small (close to the grid) scales, since the use of specialized 
front tracking algorithms'*, or correcting the numerically mea- 
sured growth rates for the observed density dilution produced 
by small scale mixing, or comparison to experiments which use 
miscible fluids, all give better agreement. 

It is important to validate numerical algorithms in experi- 
mentally accessible parameter regimes (such as the hydrody- 
namic RTI) so that these methods can be used with confidence 
to explore new physics in regimes which may be hard to realize 
in the laboratory. For example, there are a number of applica- 
tions where magnetic fields may play an important role in the 
linear evolution and nonlinear saturation of the RTI. The ax- 
ial compression of plasma produced by the ablation of wires 
in Z-pinch experiments^ is subject to the magnetic RTI. Since 
the instability can limit the maximum compression achieved in 
the pinch, understanding and controlling it is of critical impor- 
tance. Furthermore, since most astrophysical plasmas are mag- 
netized, the RTI associated with accretion onto compact ob- 
jects supernova remnants^, magnetic flux emerging from the 
solar photosphere^, and at the surface of synchrotron nebulae 
expanding into supernova ejecta^ is strongly influenced by the 
presence of magnetic fields. 

For the ideal case introduced above, the Hnear growth rate n 
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of the RTI with a uniform magnetic field of strength B parallel 
to the interface is given by' 

(B • k)2 



n =gk 



Ph-pi 



(3) 



^Ph + Pl 2TTgk(ph + Pl) 

where instabiUty occurs when > 0. Here and throughout we 
have chosen a system of units in which the magnetic permeabil- 
ity /i = 1 . In the above, k is the wavevector of a perturbation, 
and A:^ = k • k. For perturbations perpendicular to the magnetic 
field the linear growth is identical to the purely hydrodynamic 
case; these modes are often referred to as interchange modes. 
For perturbations parallel to the field, wavelengths smaller than 
the critical value = 2n /kc = / (ph - pi)g are stable, and the 
growth rate of modes at all scales larger than Ac is reduced com- 
pared to the non-magnetic case. Equivalently, for instability 
to occur on scales smaller than L, the magnetic field must be 
smaller than the critical value 

B, = \{ph-pi)gL-\''l^. (4) 
Since the growth rate is zero at both large and small k, there 
must be a wavelength Amax where the growth rate is maximum; 
it occurs at Amax = 2 Ac. The presence of a critical wavelength, a 
peak in the growth rate at Amax. and the anisotropic nature of the 
dispersion relation for perturbations parallel versus perpendic- 
ular to field lines suggests the nonlinear evolution of the mag- 
netic RTI will be much different than the non-magnetic case. 

Due to the anisotropic nature of the linear modes, it is crit- 
ical to study the magnetic RTI in full three dimensions. Two- 
dimensional ideal MHD studies in which the magnetic field is 
perpendicular to the domain are in fact equivalent to 2D hydro- 
dynamics with the gas pressure P replaced by the total pressure 
P* = f + (B • B) /2. On the other hand, two-dimensional studies 
in which the magnetic field is in the plane of the computation 
miss the interchange modes, which artificially suppresses the 
instability. For example, in 2D with B > Be the interface is 
completely stable, whereas in 3D it will be strongly unstable 
due to the interchange modes (which will act hke 2D hydrody- 
namic RTI in the plane perpendicular to the field). 

There have only been a few previous investigations of the 
magnetic RTI in full three dimensions. Jun et al.'*^ presented a 
few 3D simulations, although at much lower resolution than is 
currently possible. Their focus was in potential field amplifi- 
cation and dynamo action produced by turbulent mixing in the 
saturated state, thus most of their 3D simulations began with a 
field weak compared to Be so that Ac was initially unresolved. 
Zhu et al.'^ performed simulations of single mode interchange 
instabilities in three dimensions at very low /? = P/Pm, where 
P and Pm are the gas and magnetic pressures respectively, with 
the goal of testing the predictions of perturbation theory on the 
nonlinear structure that emerges for a single mode. More re- 
cently, Isobe et al.*^ showed that the magnetic RTI can produce 
filamentary structures during the buoyant rise of flux in the so- 
lar photosphere. 

In this paper, we use numerical MHD simulations, using 
methods that have been validated in the hydrodynamic regime, 
to study the nonlinear evolution of the magnetic RTI with strong 
fields in three dimensions. By strong we mean Ac is comparable 



to the size of the computational domain. Our goal is to explore 
how strong fields affect the formation, structure, and evolution 
of bubbles and fingers in the nonlinear regime, and how they 
affect the turbulent mixing layer (and vice-versa). Since we do 
not use front tracking methods, our calculations are hmited by 
the numerical diffusion that occurs near the grid scale^'*. By 
comparison of hydrodynamic and MHD simulations computed 
with the same parameters, numerical resolution, and algorithm, 
we can assess the relative rate of mixing between the magnetic 
and non-magnetic RTI. Interestingly, we find that even with an 
initial field too weak to resolve Ac (so that one might expect 
it to evolve more hke a hydrodynamic model), the mixing rate 
between the light and heavy fluids is substantially reduced, and 
the rate of rise of the bubbles as measured by the a parame- 
ter is substantially increased with MHD. Although such fields 
may be too weak to stabilize resolved modes, they still add a 
significant "surface tension" at very small scales, which sup- 
ports the idea that the discrepancy between the experimentally 
measured value of a in the hydrodynamic RTI and numerical 
simulations that lack front tracking is due to small scale mix- 
ing. Most of our calculations are motivated by the astrophysical 
applications of the magnetic RTI, thus we study inviscid fluids 
in the ideal MHD approximation (without any explicit resistiv- 
ity) and in a planar geometry. The latter requires the thickness 
of the mixing zone between the heavy and the hght fluids be 
much smaller than the radius of curvature of the interface R, 
and that Xc/R<^ 1. The magnetic RTI associated with Z-pinch 
experiments occurs at a very low magnetic Reynolds number, 
and in a cylindrical geometry. The inclusion of non-ideal MHD 
effects'^, and the appropriate geometry, will be important for 
future studies with application to Z-pinch'^. 

The organization of this paper is as follows. In the next sec- 
tion, we describe the equations we solve, our numerical algo- 
rithm, and the initial conditions. In section 3 we present most 
of our results, while in section 4 we summarize and discuss the 
apphcation of our results. 

2. METHOD 

We solve the equations of ideal MHD with a constant vertical 
acceleration g = (0, 0, -g) 

^ + V-(pv) = (5) 



dpy 
'df 



+ \7.(pv\-'BB) + VP* = pg 



dB 

'dt 



-hVx(vxB) = 



dE 

— + V • ((£+P*)v-B(B • V)) = pv • g 

at 



(6) 



(7) 



(8) 



In these equations, P* is the total pressure (gas plus magnetic), 
and E is the total energy density, which is related to the internal 
energy density e via 

£ = e + /9(v-v)/2 + (B-B)/2. (9) 
We use an ideal gas equation of state for which P = (j- l)e, 
where 7 is the ratio of specific heats. We take 7 = 5/3. Al- 
though we are solving the equations of compressible gas dy- 
namics, we choose a sound speed which is so large that the 



3 



resulting motions are highly subsonic and nearly incompress- 
ible. Thus, varying the adiabatic index should have httle effect 
on the results reported here. 

The calculations are performed in a three-dimensional do- 
main of size LxLx 2L, where L = 0. 1 . The vertical coordinate 
spans -0.1 < z < 0.1. The interface between heavy and Ught 
fluids is initially at z = 0. The upper half of the domain (z > 0) 
is filled with heavy fluid of density Ph = 'i, while in the lower 
half (z < 0) flie density of the light fluid is = 1. The At- 
wood number is A = 1/2. The profile of the pressure is given by 
the condition of hydrostatic equilibrium, while the amplitude is 
chosen so that the sound speed in the light fluid is unity at the 
interface, thus 

P*(z) = ^-gpz+By2 (10) 

We choose ^ = 0.1. The sound crossing time in the light fluid 
at the interface is 0.1. We will report evolutionary times in 
our computations in terms of t^. Periodic boundary conditions 
are used in the transverse (x- and >'-) directions, and reflecting 
boundary conditions are used at the top and bottom. 

The magnetic field is initialized to be uniform and along 
the X-axis B = (Bo,0,0). From the linear analysis (eq. [4]), 
if So > 5c = 0.14, then there will be no unstable wavelengths 
shorter than the size of the computational domain L. We study 
the nonlinear evolution for a variety of field strengths between 
0.1 and 0.6Bc. For the strongest field Bq = 0.6Bc, the ratio of the 
Alfven speed to the sound speed in the light fluid at the inter- 
face is 0.024, which corresponds to /? = 1.2(c^/yA)^ « 2 x 10^. 
For die weakest field Bq = O.lBc, /? = 7.5 x Thus, we study 
a regime where the magnetic energy density is small compared 
to thermal pressure (high /3), and the vertical hydrostatic equi- 
Ubrium is determined by gas pressure alone. Nonetheless, we 
study strong fields in the sense that modes paraUel to B are 
nearly completely suppressed. The magnetic RTI in plasmas 
with /3 « 1 , which is often referred to as the Parker instability in 
astrophysics'"*, has been extensively studied in the hterature'^. 

An important dimensionless parameter which characterized 
our simulations is the ratio of the critical unstable wavelength 
to the size of the computational domain Xc/L. For the strongest 
fields studied here, Xc/L = 0.36. Thus, only a few unstable 
modes are present in the domain in this case. This regime is 
appropriate to physical systems in which some scale in the sys- 
tem (e.g. the diameter of a wire in a Z-pinch) is comparable to 
the critical unstable wavelength. In §4 we discuss the interpre- 
tation of our results in terms of this parameter. 

The RTI is seeded by small amplitude, random, zone-to-zone 
perturbations to the vertical velocity added throughout the 
volume. The amplitude of the perturbations is smoothed to- 
ward the vertical boundaries; Vj.(z) = A()R(l+cos2Trz/L) where 
Ao = 0.005, and Ris a random number between -1 and 1. The 
peak perturbations are therefore 1% of c^. We have found 
that perturbing the vertical velocity is superior to perturbing 
the position of the interface, since the latter requires smooth- 
ing at the grid scale when the perturbation amphtude is smaller 
than a grid zone. Dimonte et al? were careful to introduce 
a spectrum of multimode perturbations which is truncated at 



high wavenumbers equivalent to 32 gridpoints, so that all linear 
modes are initiaUy resolved. Instead, our perturbation spectrum 
is white noise down to the grid scale. This may mean that at 
very early times, when modes are non-interacting, and the high- 
k modes dominate, there may be differences between our hy- 
drodynamical simulations. However, since we follow the mul- 
timode evolution deep into the nonlinear regime, where mode 
interactions should erase memory of the initial conditions, we 
do not expect this to limit comparisons at late times. 

The computations presented in this paper are computed using 
a recently developed Godunov method for compressible MHD 
that combines the piecewise parabohc method'^ and the di- 
rectionaUy unsplit comer transport upwind (CTU) integrator'^ 
with the constrained transport'^ algorithm for enforcing the di- 
vergence free constraint. A complete description of the algo- 
rithm, including the results of an extensive series of test prob- 
lems, is given in Ref. [19]. Adding source terms (vertical grav- 
ity) to a Godunov scheme requires particular care, by adding 
them to both the PPM reconstruction step, as well as to the 
transverse flux differences in the CTU integrator, we find our 
method holds the vertical equilibrium state (in which the pres- 
sure gradient balances gravity) exactly. 

When run in hydrodynamic mode (which utilizes the same 
algorithms as when the code is configured for MHD, except 
for the Riemann solver), the algorithm is similar to the FLASH 
and WP/PPM codes used by Dimonte et al, except for the use 
of the unspht CTU integrator. All of the simulations use a grid 
of 256 X 256 x 512, the highest resolution reported in Dimonte 
et al. In fact, apart from the perturbation spectrum, we use the 
same grid and parameters in order to more easily facihtate com- 
parisons. We have tested for convergence of our numerical so- 
lutions by running a few simulations at one half this resolution 
(128 x 128 X 256). Although such details as the number and 
location of bubbles and fingers in the flow are substantially dif- 
ferent at lower resolution, none of the diagnostics which are the 
focus of this work are changed to a significant degree. This in- 
dicates our solutions are converged with respect to these quan- 
tities. We provide a much more comprehensive investigation of 
the effect of mass diffusion at the grid scale due to numerical 
effects on our results in §3. 1 . 

The use of a fully compressible code to study low Mach 
number flows such as investigated here is not optimal, and re- 
quires many more timesteps to be taken compared to the strictly 
incompressible case (typically about 4x10"* timesteps are re- 
quired to reach f/f, = 60). The numerical algorithms used here 
are second order in both space and time^°, which helps to re- 
duce the accumulation of error. The convergence study pre- 
sented in § 3 . 1 provides insight into the effect of temporal errors 
on our results. Although approximate methods for low Mach 
number flows might be more cost effective, it is not clear they 
wiU be substantially more accurate. 

3. RESULTS 

We describe the results from four MHD simulations using 
field strengths of Bq = 0.1, 0.2, 0.4, and 0.6 B^. We shall re- 
fer to each of these calculations as runs Rl, R2, R4, and R6 
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respectively. For comparison purposes, we also describe the re- 
sults of a hydrodynamic calculation computed with the same 
code, hereafter referred to as run RH. Table 1 lists important 
properties of the calculations. In particular, note that for the 
weakest field simulation, Rl, the critical wavelength is only 
about 2.5 grid cells, therefore it is unresolved. For the smallest 
well resolved modes (requiring about 16 grid cells per wave- 
length), the growth rate is reduced by only about 15% at this 
field strength compared to the non-magnetic case. Thus, we 
expect run Rl to be similar to the purely hydrodynamic cal- 
culation RH. Our discussion will focus on the strong field run 
R6, the weak field run R2 (the weakest field for which the crit- 
ical wavelength Ac is resolved), and the hydrodynamic run RH. 
Each calculation is continued until the rising bubbles or sinking 
fingers reach the vertical boundaries, which typically requires 
about 60ts. 

3.1. Convergence Study in Two Dimensions 

Before presenting the complex nonlinear evolution of the 
magnetic RTI in fully three dimensions, it is important to first 
assess the extent to which mixing of mass and momentum due 
to numerical effects at the grid scale affect our results. We have 
used a series of simulations of the growth of the magnetic RTI 
in two dimensions, computed with the identical algorithms and 
using the same parameters as used for the three-dimensional 
runs listed in Table I, to investigate the convergence of our re- 
sults. The calculations are performed in the x-z plane, so that 
the magnetic field is in the plane of the computation. We per- 
turb the interface with a single mode with a wavelength equal to 
the horizontal size of the domain L, to investigate how numeri- 
cal grid effects alter the evolution of a single, smooth interface. 

To track the amount of mixing between the heavy and Ught 
fluids, we define a mixing parameter 

e=4/,/, (11) 

where fu is the fraction of each grid cell occupied by the heavy 
fluid (of density /?/,), and // = 1 - //,. For a purely incompress- 
ible flow, with p/, = 3 and pi = 1, then fh = (p- l)/2. In regions 
with no mixing, = 0, while the maximum value occurs for 
ifh) = ifi) = 1/2^ when 9 = 1. A useful diagnostic is the vol- 
ume averaged mixing parameter, which in two dimensions is 

{e)v = Jj4fhfidxdz/2L' (12) 

A comparison of the time evolution of (6)y measured in our 
simulations with simple analytic expectations allows us to mea- 
sure the effect of mixing at the grid scale. 

Figure 1 presents images of the mixing parameter in the evo- 
lution of the single mode RTI in two dimensions using a strong 
field (run R6) at time f/f, = 30 for four different resolutions 
corresponding to 32, 64, 128, and 256 grid points per L. Note 
the highest resolution is the same as is used in all of the three- 
dimensional results presented in this paper The images show 
that is non-zero only at the interface, which is smooth and 
has the same shape at every resolution. Clearly, the growth rate 
and nonhnear structure of a single mode is captured correctly 
even at the lowest resolution. Although the physical width of 



the interface over which mass mixing occurs is larger at low 
resolution, this is simply due to the increase in the size of grid 
cells 5x = 5z; it has not affected the rate of growth or structure of 
the mode. Figure 1 gives the visual impression that the amount 
of mass diffusion at the grid scale converges at first order (hnear 
in Sx), we quantify this dependence below. 

The initial conditions used in all of our simulations consist 
of a density (and temperature) discontinuity at the interface be- 
tween the heavy and light fluids. Initially this discontinuity is 
aligned with the grid. For the numerical algorithms used in this 
paper (Godunov method with a Roe solver), there is no numer- 
ical diffusion of this discontinuity as long as it remains at rest 
parallel to the grid. (We have confirmed that if the interface is 
not perturbed, our code holds the initial equilibrium state, and 
the mixing parameter remains Q = everywhere.) Once the in- 
terface is perturbed and the RTI begins to grow, however, mo- 
tion of the interface produces mass diffusion at the grid scale. 
This is an inevitable consequence of using a control volume ap- 
proach without interface tracking: when the interface crosses 
the middle of the cell, the volume averaged density contains 
contributions from both the heavy and light fluids, and will 
therefore be intermediate between the two. Figure 1 demon- 
strates that for the algorithms used here, the spread of the con- 
tact discontinuity is confined to a few cells (several Sx), Hence, 
at early times (when the length of the interface is just Lj,) we 
expect that 



{e)y. 



LxM5z~_MQ 



(13) 



where M is the number of cells over which the interface spreads 
due to numerical effcts, O is the normalized average of O over 
the width of the mixing region, and Ny is the number of grid 
cells in Ly . If the variation of the heavy and light fluid fractions 
is hnear over the width of the mixing region, then = 2/3. As 
the RTI grow, the length of the interface grows (e.g. figure 1). 
Thus, equation 13 expresses the expectation that the time evolu- 
tion of (O) V should be proportional to the time evolution of the 
length of the interface, and that at any time (O)v' should con- 
verge with 5x at first order (as expected for any discontinuous 
solution using a fixed grid). 

Figure 2 plots the time evolution of (0)v for the single mode 
RTI in two dimensions (run R6) at the same resolutions shown 
in figure 1. There is a rapid rise in (0)v to f/f, = 5, as the in- 
terface begins to move across the grid, and numerical diffusion 
causes it to spread to a size of a few 5x. The amplitude of (0)v 
at t/ts = 5 is in excellent quantitative agreement with the expec- 
tation of equation 1 3, if the width of the interface M « 2. There- 
after, (0)v shows slow hnear growth as the interface lengthens. 
At t/ts = 20, the mode begins to go nonhnear, the length of the 
interface begins to grow more rapidly, and (0)v increases more 
rapidly. The fractional rate of increase il/{<d)v)d{<d)v/dt is 
the same in each case, indicating the increase in time is simply 
due to the lengthening of the interface. At t /f j = 20, the con- 
vergence rate of (0)v is 0.8, close to our expectation of first 
order. 

The analysis presented thus far has considered an interface 
which remains smooth, unaffected by secondary KH instabih- 
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ties which are present with weaker magnetic fields. By twisting 
the interface on small scales, these secondary instabilities can 
increase the mass mixing at the grid scale. Figure 3 presents im- 
ages of the mixing parameter 6 at f /f , = 30 in the evolution of a 
single mode RTI computed with a resolution of 256 grid points 
per L, but with a variety of field strengths corresponding to runs 
R6 (strong field), R2 (intermediate field), Rl (weak field), and 
RH (hydrodynamic); see Table 1. The growth of an increas- 
ingly larger number of vortices at the interface is evident as the 
magnetic field strength is decreased. As these vortices wind 
up the interface, mass mixing occurs due to numerical effects 
when the radius of curvature approaches the grid scale. 

Figure 4 plots the time evolution of the volume averaged 
mixing parameter (6)v for these runs. Interestingly, until the 
appearance of the first KH roll at f / f , = 20, the amplitude and 
evolution of is identical. Thereafter, the values of 
increase rapidly in runs R2, Rl, and RH, reflecting the in- 
creased importance of the smaU scale distortion of the interface 
due to KH instabilities. 

The fact that the time evolution of (6)v is identical regard- 
less of the field strength for t j f , < 20 (until the first KH roll 
forms) indicates the intrinsic spread of a smooth interface in our 
numerical methods is independent of the field strength. This 
gives us confidence that the relative amount of mixing we ob- 
serve in three dimensions between hydrodynamic and MHD 
simulations with different field strengths will be due to phys- 
ical differences in the amount of geometric distortion of the 
interface, rather than a change in the intrinsic numerical mixing 
inherent in the algorithm between hydrodynamics and MHD. 
Moreover, increasing the resolution will not reduce this mix- 
ing, rather it will simply introduce more small scale structure 
that produces a similar amount of mixing, unless a physical pro- 
cess such as surface tension or viscosity is introduced to create 
a scale on which such motions are suppressed. We conclude the 
relative rate of mass mixing we observe between hydrodynamic 
and MHD is robust. 

The analysis above has focused on quantifying the amount of 
mass mixing due to numerical effects in our algorithms. Since 
we do not include an explicit viscosity, momentum diffusion 
at the grid scale is also dominated by grid effects. For the al- 
gorithms used here, it has been shown^° that the properties of 
the numerical diffusion of momentum produces proper conver- 
gence to a solution of the Navier-Stokes equations. The ef- 
fective Reynolds number of the flow is determined by the grid 
resolution. Since we adopt the same resolution for all runs, 
the simulations reported here are equivalent to a study of the 
change in the flow due to the effect of the magnetic field at 
fixed Reynolds number. Comparison to a specific experiment^\ 
rather than comparative results as reported here, would require 
simulations that achieve the same Reynolds number as the ex- 
periment. 

Finally, the analysis presented above has focused on the evo- 
lution of the RTI in two dimensions. In fully three dimensions, 
the value of the volume averaged mixing parameter (0)v will 
be given by equation 1 3 with the length of the interface replaced 
by the surface area. The time evolution of (O) ywiU then depend 



on the rate of growth of this surface area, which wiU be more 
rapid than in two dimensions. 

3 .2. Three-dimensional structure 

Figure 5 is a comparison of the three-dimensional structure at 
both early {t/ts = 29.6) and late {t/ts = 60) times from runs RH, 
R2, and R6 using vertical slices of the density at the edges of the 
computational domain, as well as a horizontal slice at the mid- 
plane z = 0. The hydrodynamic calculation RH (top two panels) 
can be compared directly to the results in Ref. |3|. The classic 
evolution of the hydrodynamic RTI is evident; at early times 
most structure is at small scales. Bubbles of light fluid have de- 
tached from the interface, and are dominated by a "mushroom 
cap" appearance. At late times, the dominant structures are on 
larger scales. Secondary Kelvin-Helmholtz instabilities have 
produced vortices and mixing along the edges of the bubbles, 
and a fully developed turbulent mixing layer is evident. 

In the weak field calculation R2 (middle panels), the crit- 
ical wavelength at which the magnetic field can suppress the 
RTI is small, only O.OIL, much smaller than the largest bubbles 
seen at late times in RH. Thus, we might expect the structure 
produced in the nonlinear regime in R2 to be similar to the hy- 
drodynamic case RH. At early times, R2 does show the classic 
bubble and finger morphology of the hydrodynamic RTI. There 
is no evidence for anisotropy; structures are identical parallel 
and perpendicular to the field. However, even at the early time 
there is clearly much less mixing between the light and heavy 
fluids; the horizontal slice shows that at the midplane most of 
the fluid is either at the highest or lowest density, whereas in RH 
most of the fluid is at intermediate (grey) densities (for online 
version in color, high density is red, low density is blue, inter- 
mediate density is green). At late times, the structure of the 
RTI is much different in R2 in comparison to RH. Rather than 
a turbulent mixing layer, in R2 the bubbles have become long 
columns which have extended far above the midplane. These 
structures show no anisotropy. There continues to be little mix- 
ing between the fluids. Thus, even a weak field has strongly 
affected the structure. 

In the strong field calculation R6 (bottom panels) the 
anisotropic nature of linear modes is clearly evident: the density 
fluctuations at the midplane are aligned with the direction of the 
field (along the x-axis). As in R2, very little mixing is evident 
at the midplane in R6. At late times in the strong field case, 
the dominant structures are smooth and highly anisotropic. The 
slice in the y-z plane shows colunms and bubbles which result 
from interchange modes that act like the hydrodynamic RTI in 
two dimensions. In the x-z plane, however, extended loop-like 
structures are evident which result from suppression of short 
wavelength modes along the field. These structures are unlike 
any of the bubbles seen in the hydrodynamic RTI (run RH, top 
panel). 

To further illustrate the nonlinear structure at late time. Fig- 
ure 6 shows isosurfaces of the density at p = 1.1 and 2.9, along 
with vertical slices of the density at the faces of the computation 
domain for runs RH and R2 at f//, = 56. A complex network 
of bubbles is evident in the shape of the density isosurface in 
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RH, whereas in R2 the bubbles are much larger and smoother. 
Note the circular ring near the rear edge of the domain in RH 
resulting from the roll-up of the KH instability around a spheri- 
cal bubble. The slices at the edge of the domain also make clear 
the turbulent mixing in RH, whereas in R2 there is far less mix. 

Isosurfaces and slices of the density in the strong field run R6 
are shown at two times in Figure 7; the left panel is at f/fj = 20 
and the right at / //j = 56. At the earher time, the isosurface 
reveals the formation of filaments and tubes rather than bub- 
bles and fingers as in the hydrodynamic case. The anisotropy 
of linear modes is evident in the comparison of the slices in the 
x-z and >'-z planes. In the former, the magnetic field has sup- 
pressed short wavelengths, and the density structures are highly 
elongated in the x— direction. In the latter, the magnetic field 
is perpendicular to the plane, and thus has no effect. Bubbles 
and fingers on short wavelengths have formed, reminiscent of 
the 2D hydrodynamic RTI. The formation of flux tubes seen 
at early times is very similar to that observed in the solar pho- 
tosphere, as shown in Ref. [8]. At late times (right panel of 
Figure 7), the density isosurface is remarkably smooth. Rather 
than sheets, more rounded bubbles have been formed at the tips 
of the columns by the flow of fluid along loops until it collects 
at the tips. Similar evolution is observed in the nonlinear regime 
of the Parker instability'^. 

An important diagnostic of the instability is the rate at which 
bubbles and fingers are displaced from the initial interface. To 
quantify this rate, we define the horizontally averaged mass 
fraction of heavy fluid as (following Ref. [3]) 

(/,) = / / fhdxdy/L\ (14) 

JxJy 

Since our simulations are not purely incompressible, the maxi- 
mum (minimum) densities can be slightly larger (smaller) than 
3 (one). To account for this, we define the height of bubbles to 
be the location where {fh) = 0.985, while the height of fingers 
is the location where (//,) = 0.015. 

Figure 8 is a plot of the height of bubbles and fingers in runs 
RH, R2, and R6 versus t^. From the self-similar arguments that 
lead to eq. [1], we expect at late time a straight line with a slope 
of a. In each case, this expectation is confirmed. In the hydro- 
dynamic run RH, the best fit slope at late time is a = 0.021, 
whereas for MHD we obtain a = 0.035 for both R4 and R6. Al- 
though not plotted to avoid cluttering the figure, we have also 
confirmed that for the weakest field run Rl, the height of bub- 
bles follows a straight line when plotted versus f', with a slope 
of a = 0.031. Thus, we find that (1) the rate at which bub- 
bles rise in the hydrodynamic RTI measured in our simulations 
agrees with the results of Dimonte et al for non-front tracking 
methods, and (2) the addition of even a small magnetic field sig- 
nificantly increases the slope. We show below that there is far 
less mixing in the MHD RTI simulations, which may account 
for this increase. 

3.3. Mixing 

Figure 9 plots the vertical profile of the fraction of heavy fluid 
{fh) for runs RH and R6 at t/ts = 56, where the horizontal axis 



has been scaled by the bubble height at that time. The profiles 
are remarkable similar to each other. We also find the profiles 
are the same at different times in the evolution, as is expected 
due to the self-similar evolution. The similarities in the profiles 
of {fh) between the magnetized and unmagnetized RTI, despite 
the lack of a turbulent mixing layer in the former, suggests that 
this profile is not sensitive to mixing. 

Figure 10 plots the vertical profile of the horizontally aver- 
aged mixing parameter 

(0)=4(/„//) (15) 
for runs RH, Rl, R2 and R6 at f/f, = 56. There is a monotonic 
decrease in the maximum value of (9) with the field strength. 
The hydrodynamic case RH is closest to reaching (O) = 1, the 
maximum value possible. Even the weakest field case, run Rl, 
in which the critical wavelength suppressed by the magnetic 
field is unresolved initially, has a significantly lower value of 
{&), indicating far less mixing is occurring compared to the 
hydrodynamic case. 

It also is of interest to compare the integral of the mixing pa- 
rameter over the vertical coordinate between different runs. For 
the hydrodynamic and weak field runs RH and Rl respectively, 
/ {Q)dz ~ 0.06, with smaller values obtained for the stronger 
field cases R2 and R6. This quantity has dimensions of length 
and so can be interpreted as the effective width of a completely 
mixed region (for which {Affh) = 1). The volume averaged 
value (9)v = / {0)dz/L « 0.6, which can be interpreted as 
the fractional height of the domain over which mixing occurs. 
Thus, not only is the mixing local to the initial interface larger 
for hydrodynamics compared to MHD, but also even in an in- 
tegral sense the effective width is reduced. 

We have found the variance of the density is also a sensitive 
diagnostic of the mixing region. At each vertical position z, 
we compute the horizontally averaged density {p), where the () 
denotes an average over a horizontal plane as in eq. 11. The 
variance is then 

Spiz)={ip-{p))Y^ (16) 

Regions which are well mixed have a smaller variance. If the 
heavy and light fluids remain completely unmixed (that is, if 
the density at every location can only be either 3 or one) then 
the largest value of the variance is 5p/ {p) = 1/ \/3 « 0.57 and 
occurs for (//,) = 1 /4. 

Figure 1 1 plots the vertical profile of dp/ {p) in Runs RH, R2, 
and R6 at r/ = 56. Both the weak and strong field runs are sim- 
ilar to each other, and have a much larger variance than the un- 
magnetized case, indicating less mixing. For (//,) = 1 /2, which 
from a visual inspections of the horizontal slice plane in Figure 
5 is a rough approximation to the volume fraction of the heavy 
fluid there, the ideal case of no mixing gives Sp/ {p) = 1 /2. 
From Figure 1 1 the peak values in runs R2 and R6 approach 
this value, an indication that very little mix occurs in the MHD 
RTI. Note from Figure 9 fliat {fh) = 1 /4 occurs at Z/H of -0.4. 
At the time of this plot, the height of the bubbles H « 0.7, thus 
we should expect the peak in the variance to occur at z « -0.03, 
which is in good agreement with Figure 1 1 . Thus, the asymme- 
try in the variance (with larger values occurring for negative z) 
is a consequence of the lower (//,) there. 
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To investigate the dynamical processes that can lead to re- 
duced mixing in the MHD RTI compared to hydrodynamics, 
we have calculated the distribution of the amplitude of the vor- 
ticity I W 1=1 V X v |. Figure 12 plots the number of cells A' in 
runs RH and Rl which have a given value of the amplitude of 
vorticity, A'^(| W |), versus \W \ . Clearly there is a shift in the 
distribution to lower values of | W | in the weak field run Rl 
(dotted line). Large values of the vorticity are associated with 
velocity shear over small scales. The shift in the distribution 
to lower I W \ with weak fields is an indication that even if the 
critical wavelength Ac is not resolved in the initial conditions 
(Rl), the tension forces associated with small scale bending of 
the field lines is sufficient to affect the flow, and suppress the 
shear. We postulate that the suppression of small scale shearing 
motion implied by the change in the distribution A^(| W |) is re- 
sponsible for the reduced mixing observed in the MHD RTI as 
shown by figures 10 and 1 1 . 

3.4. Evolutionof the magnetic field 

To investigate the three-dimensional structure of the mag- 
netic field, Figure 1 3 plots sUces of the magnetic energy in the 
fluctuating part of the field, 

^BV2 = (B,-Bo)V2+Bv/2+B^/2 (17) 
at the edges of the domain, and at the midplane z = in runs R2 
and R6 at f/fj = 60. The slices are made at the same locations 
and the same times as the sUces of the density shown in the 
right hand panels of figure 5. There is a direct correspondence 
to the structures visible in the magnetic field and the density. 
The fluctuations in 8B^ occur at smaller scales in R2 in com- 
parison to R6. Most of the magnetic energy is concentrated in 
rope and sheet-like structures associated with the bubbles and 
fingers in the density. The maximum of the magnetic energy 
is quite similar in the two plots (the maximum is 0.021 in R2 
and 0.015 in R6), despite the fact the energy in the background 
field is nearly an order of magnitude larger in R6. In fact, at 
t/ts = 60, the energy in fluctuations is larger in the weaker field 
run. As discussed below, this is Ukely due to the fact that the 
fingers and bubbles are larger in R2 at this time, thus more grav- 
itational energy has been released that can be tapped to amplify 
the magnetic field. 

Figure 14 is a plot of the the vertical profile of the horizon- 
tally averaged magnetic energy at two times in both runs R2 and 
R6. The energy in the vertical component {B\/2) is separated 
from the horizontal components {{B^-Bof- /2+By/2). The en- 
ergy in the vertical field always dominates, a consequence of 
the dominance of vertical motions. The amplitude of the en- 
ergy in the vertical field is remarkably similar between the two 
runs; ait/ts = 56 it is about 0.009 in R2 and 0.0075 in R6. This 
is consistent with the amplitude of the energies being similar in 
the slices shown in figure 13. The total energy in fluctuations 
(the integral over vertical height of the Unes plotted in figure 
14) must therefore be similar between the two runs; we discuss 
this further below. The ratio of energy in the vertical versus 
horizontal fields is larger in the weak field case R2. 

Simple energy arguments can be used to interpret the am- 
phtude of the magnetic energies observed in figures 13 and 



14. The gravitational potential energy released by descending 
heavy fluid is converted into kinetic energy by the RTI and sec- 
ondary KH instabiUties. In turn, these result in twisting and 
amplification of the magnetic field. Some kinetic and magnetic 
energy is converted into internal energy by viscous and resistive 
dissipation. Although our simulations do not include expUcit 
dissipation, our method conserves total energy exactly so that 
whatever energy is lost by numerical viscosity and reconnec- 
tion is expUcitly captured as an increase in internal energy. The 
simplest expectation is that the kinetic and magnetic energies 
will remain in approximate equipartition, and that the amount 
of energy released when the tip of the fingers reaches the same 
height (as measured by the vertical location where (/;,) = 0.995) 
will be the same amongst all runs. To test these ideas, table 2 
compares the volume averaged energies from all the runs when 
the tips of the fingers reaches z/H = 0.5. The third column gives 
the time at which this happens, the fourth is the magnetic en- 
ergy in fluctuations, the fifth is the kinetic energy, the sixth is 
the change in internal energy at that time. The seventh column 
is the total change in energy, and the last is the amplification 
factor of the magnetic energy in fluctuations. Since the data is 
collected not at the same physical time, but at the point where 
the fingers have reached the same height (and therefore the 
same amount of gravitational energy has been released), then 
the total of the kinetic, magnetic, and thermal energies should 
be the same amongst all the runs. Indeed, this expectation is 
confirmed, the difference in total energy is only 15% amongst 
all the simulations. Note also the magnetic and kinetic energies 
are in rough equipartition, regardless of the initial field strength. 
The increase in internal energy is a non-negligible contribution 
to the total change, especially in the weaker field simulations. 
In fact the total gravitational potential energy released depends 
on the detailed distribution of density in the interface region 
and not just the location of the tip of the fingers, thus we do 
not expect the total energy to be identical between the different 
runs. 

These same energy arguments can be used to predict the time 
evolution of the magnetic and kinetic energies. The gravita- 
tional energy released depends on the total mass and the dis- 
tance it drops, both of which are proportional to the height of 
the bubbles h. Thus, the energy released E oc (x t'^. Figure 
15 plots each component of the magnetic and kinetic energies 
in runs Rl, R2, R4, and R6 versus t'^. Our expectation is that 
each should be a straight line, with rough equipartition between 
the transverse magnetic and kinetic energies, and with the ver- 
tical components dominating. This expectation is clearly borne 
out by the figure, at least so long as the energy in the vertical 
component of the field B^/2 (which grows the fastest) is less 
than the magnetic energy associated with the initial field Bq/2. 
Thus, during the evolution of strong field case R6 (figure 15a), 
the ratio B^ /By is always less than one, and the growth of each 
component of the energy closely follows t'^. With increasingly 
weaker fields, the amplification of B^/Bq becomes larger and 
larger, and the time evolution of each component diverges far- 
ther and farther from z'*. For the very weak field case, Rl (figure 
15d), the ratio B^/Bq = 1 is reached very early in the evolution. 
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and the maximum amplification at late time is more than 15. 
In this case, the time evolution of the energy is lower than t'^. 
As discussed in the next section, we expect deviation from the 
simple f"* scaling when the bubbles and fingers have moved a 
distance much larger then Ac. 

4. SUMMARY AND DISCUSSION 

We have studied the RTI in ideal MHD and full three dimen- 
sions. We have restricted our analysis to uniform fields parallel 
to the interface. We study the high-/3 regime, where the energy 
density in the magnetic field is small compared to the thermal 
energy in the fluid. Nonetheless, we study strong fields in the 
sense that the initial field strength Bq < Be, where Be is the 
critical field strength at which all modes parallel to B are com- 
pletely suppressed. We use numerical methods that have been 
validated in the sense that they reproduce the growth rate of fin- 
gers and the amount of mixing between light and heavy fluids in 
the hydrodynamic RTI, as reported in previous high resolution 
numerical simulations^, and laboratory experiment^^^. 

Uniform magnetic fields do not suppress the RTI, but rather 
make modes strongly anisotropic. Along the field, the growth 
rate of modes is reduced, and wavelengths below a critical value 
are stable. Perpendicular to the field, the dispersion relation 
of the magnetic RTI is identical to the dispersion relation in 
hydrodynamics. Even if the field is made arbitrarily strong to 
suppress all modes parallel to B, the density interface will still 
be RTI unstable in 3D due to the interchange modes perpendic- 
ular to B. It is therefore critical to study the MHD RTI in full 
three dimensions. 

Although uniform magnetic fields can not suppress the RTI 
in three dimensions, they significantly change the nonlinear 
evolution and structure. Even weak fields reduce the mix be- 
tween light and heavy fluids, resulting in fingers and bubbles 
which rise much more quickly compared to the hydrodynamics 
case. Such fields are weak in the sense that Xc/L is small, how- 
ever they still are strong enough to influence the flow through 
tension forces at small scales, as evidenced by the change in 
distribution of the vorticity between very weak field, and purely 
hydrodynamical, simulations (figure 12). A turbulence mix- 
ing zone is not produced with strong fields. Instead, at early 
times the bubbles and fingers are elongated along B, forming 
flux ropes and tubes. Fluid drains along these tubes, pooling 
as bubbles at the tips and vaUeys, eventually forming the usual 
bubble and finger morphology. Interchange instabiUties wrin- 
kle the surfaces of the bubbles perpendicular to the field. Sev- 
eral diagnostics, including the variance of the density, that are 
good diagnostics of the amount of mixing between heavy and 
Ught fluids show that there is a monotonic decrease in mixing 
with increasing field strength. 

The evolution of the MHD RTI follows the same self-similar 
evolution as in hydrodynamics. The vertical profile of the vol- 
ume fraction of heavy fluid is self-similar. The vertical extent 
of the bubbles and fingers h cx f^. Simple energy arguments 
suggest the magnetic and kinetic energies should grow in time 
as t'^, and our results confirm this expectation. The energy in 
transverse motions and field are always in rough equipartition. 



The total energy in fluctuations is independent of the initial field 
strength, but depends only on the extent of the bubbles and fin- 
gers (which in turn determines the amount of gravitational en- 
ergy released that is available for field amplification). Since 
the magnetic energy density in the mixing layer grows to the 
same value in every case (determined by the amount of gravi- 
tational binding energy released by descending heavy fluid and 
the thickness of the mixing layer), this leads to a much larger 
fractional increase in magnetic energy for initially weak fields 
(more than a factor of 10). However, this dynamo action can 
never lead to strong fields with /3 ^ 1 , since the total magnetic 
energy is Umited to equipartition with the kinetic rather than 
thermal energy, and the flows induced by the RTI are highly 
subsonic. One explanation for the larger increase in magnetic 
energy with initially weak fields is that the critical wavelength 
at which the RTI is stable is smaller for weak fields, thus weak 
fields can be folded, twisted, and amplified on smaller scales 
than strong fields. 

The self- similar evolution of the magnetic and kinetic ener- 
gies diverges from the simple expectation above once the mag- 
netic energy in the vertical component of the field exceeds that 
associated with the original field, i.e. B^/Bq > 1. This oc- 
curs when the bubbles and fingers have moved a large distance 
compared to Ac. The evolution of the energies in the strong 
field simulations described here (R6), which has Xc/L = 0.36, 
closely follows a f scaling throughout. However, for the very 
weak field simulation (Rl), which has Xc/L = 0.01, this scahng 
is broken at early times. This difference can be interpreted as 
being due to the ve ry dif ferent dimensionless length and time 
scales, Xc/L and \/Xc/g/ts respectively, in each case. If the 
very weak field case Rl were repeated with identical parame- 
ters but in a computational domain of size L/36, then it would 
have the same ratio Xc/L, and it would evolve identically to the 
strong field case R6 on a time scale which is 6 times shorter. 
Similarily, if the strong field case R6 were repeated in a com- 
putational domain 36 times larger, it would be identical to run 
Rl at times a factor of 6 longer than Rl. Thus, run Rl samples 
a much later stage of evolution of the magnetic RTI than R6. 
We conclude the self-similar evolution which predicts energy 
growth at is only applicable for t < Tc= \J Xc/g. Although 
we have reported simulations of the magnetic RTI with differ- 
ent field strengths and the same sized computational domain L 
in this paper, this is equivalent to simulations with the same 
field strength but with different sized computational domains. 

The fact that the saturated magnetic energy is independent 
of the initial field strength provides a simple explanation for 
the break in the t'^ scaling of the magnetic energy at late times 
Tc. In this case, since the magnetic field strength and mag- 
netic energy density B^ reaches a constant value in the mixing 
layer, then the total energy wiU increase further only because of 
the thickening of the mixing layer, thus in the saturated regime 
E Q(.h(xt^. 

The MHD RTI is relevant to a number of astrophysical sys- 
tems, for example density interfaces in the ISM, to the forma- 
tion of filaments in the Crab nebulae^, and to the contact discon- 
tinuity in supemovae blast waves'. As discussed above, a uni- 
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form field will not suppress the RH in these systems, nor even 
reduce the growth rate of modes perpendicular to B. It wiU, 
however, result in highly anisotropic structures, and reduce the 
mixing between fluids. As reported elsewhere (Stone & Gar- 
diner 2007, submitted to the Astrophys. J.), we have also found 
that rotation of the direction of the field with vertical position 
(which is equivalent to currents parallel to the interface) signif- 
icantly affects the nonhnear evolution of the the MHD RTI. 

The MHD RTI is also relevant to Z-pinch experiments, where 
implosion is driven by low density, highly magnetized plasma 
ablated from wires. Non-ideal MHD effects (resistivity, and 
Hall currents) are important at the densities and temperatures 



realized in the plasma in these experiments^. Studies of the 
MHD RTI including resistive dissipation and Hall currents in a 
cylindrical geometry with higher (3 are needed^^. 
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Table 1 
Parameters of Runs 



Run B/B, Ac/L AmaxA -^c/Ay 



RH 


0.0 








Rl 


0.1 


0.01 


0.02 


2.56 


R2 


0.2 


0.04 


0.08 


10.2 


R4 


0.4 


0.16 


0.32 


41.0 


R6 


0.6 


0.36 


0.72 


92.2 



Table 2 

Volume averaged energies at h/L = 0.5. 



Run bp /be tju e-ep total 5b'^/bl 



RH 




53 




1.6 


2.5 


4.1 




Rl 


0.1 


46 


0.85 


1.6 


2.7 


4.4 


10.6 


R2 


0.2 


42 


1.0 


1.3 


2.1 


4.4 


3.2 


R4 


0.4 


42 


1.2 


1.4 


2.1 


4.7 


0.93 


R6 


0.6 


49 


1.1 


1.2 


2.2 


4.5 


0.37 



11 




Fig. 1. — Images of the mixing parameter 0, defined in equation 12, in the two-dimensional version of run R6 at f/f, = 30.0 at resolutions of 32 X 64 (left), 
64 X 128 (middle left), 128 X 256 (middle right), and 256 X 512 (right). The color table runs blue to red (online version) over the range zero to one, with white 
corresponding to = 0. 
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Fig. 2. — Time evolution of tlie volume-averaged mixing parameter at different resolutions for the two-dimensional version of the strong field run R6. Each curve 
is labelled by the number of grid points per L. 
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Fig. 3. — Images of the mixing parameter © at t/ts = 30.0, in two-dimensional version of runs R6 (left, strong field), R2 (middle left, intermediate field), Rl 
(middle right, weak field), and RH (right, hydrodynamic), all at a resolution of 256 X 512. The color table runs blue to red (online version) over the range zero to 
one, with white corresponding to © = 0. 
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Fig. 4. — Time evolution of ttie volume-averaged mixing parameter for runs R6 (strong field), R2 (intermediate field), Rl (weak field) and RH (hydrodynamic). 
The rapid increase in the mixing rate occurring ait/ts = 20 in runs with weaker fields is due to the formation of KH rolls, as is evident in figure 3. 
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Fig. 5. — Slices of the density aXtjts = 29.6 (left panels) and ?/fj = 60.0 (right panels) in rans RH (top, hydrodynamic), R2 (middle, weak field), and R6 (bottom, 
strong field). Note the decrease in mixing in the MHD RTI (as evidenced by reduced volume at intermediate densities, i.e. grey regions), and the elongation of 
structures along the magnetic field (x-axis) in the strong field case. (Online version in color.) 
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Fig. 6. — Isosurfaces of the density at p = 2.9 and p = 1.1 at t jts = 56.0 in runs RH (left, hydrodynamic) and R2 (right, weak field). Also shown are slices of the 
density at the edges of the domain. (Online version in color) 
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Fig. 7. — Isosurfaces of the density at p = 2.9 and p = 1.1 at t/ts = 20.0 (left) and t/ts = 56.0 (right) in run R6 (strong field). Also shown are slices of the density 
at the edges of the domain. (Online version in color.) 
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Fig. 8. — Height of bubbles (z > 0) and fingers (z < 0) as a function of time in runs RH (solid Une, hydrodynamic), R2 (dashed line, weak field) and R6 (dotted 
line, strong field). The height is defined as the location where the horizontally averaged mass fraction (eq. 1 1) is 0.985 and 0.015 respectively. 




Fig. 9. — Vertical profile of the mass fraction (/;,) defined by eq. 11 in runs RH (dashed line, hydrodynamic) and R6 (solid line, strong field) at tjts = 56. The 
horizontal axis has been scaled by the height of the bubbles at that time. 
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Fig. 10. — Vertical profile of the horizontally averaged mixing parameter (0) = 4{fi,fi) at time t/ts =56 in runs RH (solid Une, hydrodynamic), Rl (dotted Une, 
very weak field), R2 (dot-dash Une, weak field), and R6 (dashed line, strong field). As measitfed by (0), there is monotonicaUy less mixing with increasing field 
strength. 
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Fig. 1 1 . — Vertical profile of the variance of the density defined by eq. 13 in RH (solid line, hydrodynamic), R2 (dotted line, weak field) and R6 (dashed line, 
strong field). Both of the MHD runs show a significantly larger variance than the hydrodynamic case indicating less mixing. 
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Fig. 12. — Distribution of tlie amplitude of the vorticity W in runs RH (solid line, hydrodynamic) and Rl (dotted line, very weak field) at t/ts 
axis is the number of cells with the given value of W. The horizontal axis is scaled by the maximum of W in run RH. 
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Fig. 13. — Slices of the magnetic energy in fluctuations, defined by eq. 14, at ///,, = 60 in runs R2 (left, weak field) and R6 (right, strong field). The slices are 
taken at the same locations and at the same time as the right-hand panels in figure 1. (Online version in color.) 
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Fig. 14. — Profiles of the horizontally averaged magnetic energy in fluctuations for runs R2 (top, weak field) and R6 (bottom, strong field). The dashed lines in 
each plot correspond to the energy in the vertical component of the field B^, while the solid lines correspond to the the energy in the horizontal components of the 
field (Bx -Bq)^ +By. The profiles are shown at t/tj = 28 and t/ts = 56 in each plot, the latter pair of lines extend over a larger horizontal range in both plots. 
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Fig. 15. — Time evolution of each component of the magnetic and kinetic energies in runs (a) R6, strong field, (h) R4, intermediate field, (c) R2, weak field, 
and (d) Rl, very weak field. Solid lines show magnetic energy, dashed lines are kinetic. Each line is labeled by the associated field component. The energy in 
perturbations is shown for the ;c-component of the magnetic field, SB^/l = (B^ -Bq)/2, and in each panel the values are scaled by the initial magnetic energy 



